# This code replicates Figure 10
rm(list=ls())
setwd("")
library("gridExtra")
library(foreign)
library(ri)
library(ggplot2)
data <- read.dta("data.dta")
data <- data[complete.cases(data$exp4_swayed), ]

#### Adventist
y <- data$exp4_swayed
Z <- data$treatment_adventist
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE

## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

p <- dispdist(distout, ate)[3]
ate
distout <- as.data.frame(distout)
m1 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey", bins=50) + 
      geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
      annotate("text", x = -0.06, y = 42, label = "ATE: 0.105 \ \ \ \ \ ") +
      annotate("text", x = -0.06, y = 40, label = "P-Value: 0.029") +
      labs(title = "Adventist Treatment Effect") + 
      labs(x = "Estimated ATE") +
      labs(y = "Frequency")
m1


#### Mixed

y <- data$exp4_swayed
Z <- data$treatment_mixed
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE
## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

dispdist(distout, ate)[4]
ate
distout <- as.data.frame(distout)
m2 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey", bins=50) + 
  geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
  annotate("text", x = -0.06, y = 42, label = "ATE: -0.003 \ \ \ \ ") +
  annotate("text", x = -0.06, y = 40, label = "P-Value: 0.325") +
  labs(title = "Mixed Treatment Effect") + 
  labs(x = "Estimated ATE") +
  labs(y = "Frequency")
m2


#### Peruana

y <- data$exp4_swayed
Z <- data$treatment_peruana
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE

## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

dispdist(distout, ate)[4]
ate
distout <- as.data.frame(distout)
m3 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey", bins=50) + 
  geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
  annotate("text", x = -0.06, y = 42, label = "ATE: -0.017 \ \ \ \ ") +
  annotate("text", x = -0.06, y = 40, label = "P-Value: 0.454") +
  labs(title = "Peruana Treatment Effect") + 
  labs(x = "Estimated ATE") +
  labs(y = "Frequency")
m3


#### Maranatha

y <- data$exp4_swayed
Z <- data$treatment_maranatha
cluster <- data$village
perms <- genperms(Z,clustvar=cluster)       # all possible permutations
probs <- genprobexact(Z, clustvar=cluster)   # probability of treatment
ate <- estate(y,Z,prob=probs)                               # estimate the ATE

## Conduct Sharp Null Hypothesis Test of Zero Effect for Each Unit
Ys <- genouts(y,Z,ate=0) # generate potential outcomes under sharp null of no effect
distout <- gendist(Ys,perms, prob=probs) # generate sampling dist. under sharp null
dispdist(distout, ate) # display characteristics of sampling dist. for inference

dispdist(distout, ate)[4]
ate
distout <- as.data.frame(distout)
m4 <- ggplot(distout, aes(x = distout)) + geom_histogram(colour = "black", fill = "grey", bins=50) + 
  geom_vline(xintercept = ate, colour="red", linetype = "longdash") + 
  annotate("text", x = -0.06, y = 100, label = "ATE: -0.081 \ \ \ \ \ ") +
  annotate("text", x = -0.06, y = 94, label = "P-Value: 0.007") +
  labs(title = "Maranatha Treatment Effect") + 
  labs(x = "Estimated ATE") +
  labs(y = "Frequency")
m4


library("gridExtra")
grid.arrange(m1, m2, m3, m4, ncol=2, nrow=2, 
             widths=c(3, 3), heights=c(3, 3))
pdf("persuasion_randomization.pdf", width = 9, height = 9)
grid.arrange(m1, m2, m3, m4, ncol=2, nrow=2, 
             widths=c(3, 3), heights=c(3, 3))
dev.off()

